library(tidyverse)
library(gtsummary)
library(haven)
library(skimr)
library(corrplot)
library(bkmr)
library(forestplot)
library(mice)
library(rms)
library(tableone)
library(ggplot2)
library(ggcorrplot)
getwd()
# setwd('D:/rCode/论文实战/paper3')
#获取2007至2016年所有研究对象
setwd('G:/BaiduNetdiskDownload/nhanes')
demo_e <- read_xpt("2007-2008/Demographics/demo_e.xpt")
demo_f <- read_xpt("2009-2010/Demographics/demo_f.xpt")
demo_g <- read_xpt("2011-2012/Demographics/demo_g.xpt")
demo_h <- read_xpt("2013-2014/Demographics/demo_h.xpt")
demo_i <- read_xpt("2015-2016/Demographics/demo_i.xpt")
demo_all <- dplyr::bind_rows(list(demo_e,demo_f,demo_g,demo_h,demo_i))

#获取所有实验室尿酸数据
Laboratory_e <- read_xpt("2007-2008/Laboratory/BIOPRO_e.xpt")
Laboratory_f <- read_xpt("2009-2010/Laboratory/BIOPRO_f.xpt")
Laboratory_g <- read_xpt("2011-2012/Laboratory/BIOPRO_g.xpt")
Laboratory_h <- read_xpt("2013-2014/Laboratory/BIOPRO_h.xpt")
Laboratory_i <- read_xpt("2015-2016/Laboratory/BIOPRO_i.xpt")
Laboratory_all <- dplyr::bind_rows(list(Laboratory_e,Laboratory_f,Laboratory_g,Laboratory_h,Laboratory_i))

#获取实验室中尿液邻苯二甲酸酯数据
PHTHTE_e <- read_xpt("2007-2008/Laboratory/PHTHTE_e.xpt")
PHTHTE_f <- read_xpt("2009-2010/Laboratory/PHTHTE_f.xpt")
PHTHTE_g <- read_xpt("2011-2012/Laboratory/PHTHTE_g.xpt")
PHTHTE_h <- read_xpt("2013-2014/Laboratory/PHTHTE_h.xpt")
PHTHTE_i <- read_xpt("2015-2016/Laboratory/PHTHTE_i.xpt")
PHTHTE_all <- dplyr::bind_rows(list(PHTHTE_e,PHTHTE_f,PHTHTE_g,PHTHTE_h,PHTHTE_i))

#获取BMI数据
BMI_e <- read_xpt("2007-2008/Examination/bmx_e.xpt")
BMI_f <- read_xpt("2009-2010/Examination/bmx_f.xpt")
BMI_g <- read_xpt("2011-2012/Examination/bmx_g.xpt")
BMI_h <- read_xpt("2013-2014/Examination/bmx_h.xpt")
BMI_i <- read_xpt("2015-2016/Examination/bmx_i.xpt")
BMI_all <- dplyr::bind_rows(list(BMI_e,BMI_f,BMI_g,BMI_h,BMI_i))

#获取血清可替宁数据
cotnal_e <- read_xpt("2007-2008/Laboratory/cotnal_e.xpt")
cotnal_f <- read_xpt("2009-2010/Laboratory/cotnal_f.xpt")
cotnal_g <- read_xpt("2011-2012/Laboratory/cotnal_g.xpt")
cot_h <- read_xpt("2013-2014/Laboratory/cot_h.xpt")
cot_i <- read_xpt("2015-2016/Laboratory/cot_i.xpt")
cot_all <- dplyr::bind_rows(list(cotnal_e,cotnal_f,cotnal_g,cot_h,cot_i)) 

#获取饮酒数据
alq_e <- read_xpt("2007-2008/Questionnaire/alq_e.xpt")
alq_f <- read_xpt("2009-2010/Questionnaire/alq_f.xpt")
alq_g <- read_xpt("2011-2012/Questionnaire/alq_g.xpt")
alq_h <- read_xpt("2013-2014/Questionnaire/alq_h.xpt")
alq_i <- read_xpt("2015-2016/Questionnaire/alq_i.xpt")
alq_all <- dplyr::bind_rows(list(alq_e,alq_f,alq_g,alq_h,alq_i)) 


#获取尿肌酐数据
alb_cr_e <- read_xpt("2007-2008/Laboratory/alb_cr_e.xpt")
alb_cr_f <- read_xpt("2009-2010/Laboratory/alb_cr_f.xpt")
alb_cr_g <- read_xpt("2011-2012/Laboratory/alb_cr_g.xpt")
alb_cr_h <- read_xpt("2013-2014/Laboratory/alb_cr_h.xpt")
alb_cr_i <- read_xpt("2015-2016/Laboratory/alb_cr_i.xpt")
alb_cr_all <- dplyr::bind_rows(list(alb_cr_e,alb_cr_f,alb_cr_g,alb_cr_h,alb_cr_i))

#获取血压问卷数据
bpq.e <- read_xpt("2007-2008/Questionnaire/bpq_e.xpt")
bpq.f <- read_xpt("2009-2010/Questionnaire/bpq_f.xpt")
bpq.g <- read_xpt("2011-2012/Questionnaire/bpq_g.xpt")
bpq.h <- read_xpt("2013-2014/Questionnaire/bpq_h.xpt")
bpq.i <- read_xpt("2015-2016/Questionnaire/bpq_i.xpt")
bpq_all <- dplyr::bind_rows(list(bpq.e,bpq.f,bpq.g,bpq.h,bpq.i))

#获取活动数据
paq_e <- read_xpt("2007-2008/Questionnaire/paq_e.xpt")
paq_f <- read_xpt("2009-2010/Questionnaire/paq_f.xpt")
paq_g <- read_xpt("2011-2012/Questionnaire/paq_g.xpt")
paq_h <- read_xpt("2013-2014/Questionnaire/paq_h.xpt")
paq_i <- read_xpt("2015-2016/Questionnaire/paq_i.xpt")
paq_all <- dplyr::bind_rows(list(paq_e,paq_f,paq_g,paq_h,paq_i))

#获取糖化血红蛋白数据
ghb_e <- read_xpt("2007-2008/Laboratory/ghb_e.xpt")
ghb_f <- read_xpt("2009-2010/Laboratory/ghb_f.xpt")
ghb_g <- read_xpt("2011-2012/Laboratory/ghb_g.xpt")
ghb_h <- read_xpt("2013-2014/Laboratory/ghb_h.xpt")
ghb_i <- read_xpt("2015-2016/Laboratory/ghb_i.xpt")
ghb_all <- dplyr::bind_rows(list(ghb_e,ghb_f,ghb_g,ghb_h,ghb_i))
ghb_all$LBXGH
#获取空腹血糖数据
glu_e <- read_xpt("2007-2008/Laboratory/glu_e.xpt")
glu_f <- read_xpt("2009-2010/Laboratory/glu_f.xpt")
glu_g <- read_xpt("2011-2012/Laboratory/glu_g.xpt")
glu_h <- read_xpt("2013-2014/Laboratory/glu_h.xpt")
glu_i <- read_xpt("2015-2016/Laboratory/glu_i.xpt")
glu_all <- dplyr::bind_rows(list(glu_e,glu_f,glu_g,glu_h,glu_i))
glu_all$LBXGLU
#获取OGTT数据
OGTT_e <- read_xpt("2007-2008/Laboratory/OGTT_e.xpt")
OGTT_f <- read_xpt("2009-2010/Laboratory/OGTT_f.xpt")
OGTT_g <- read_xpt("2011-2012/Laboratory/OGTT_g.xpt")
OGTT_h <- read_xpt("2013-2014/Laboratory/OGTT_h.xpt")
OGTT_i <- read_xpt("2015-2016/Laboratory/OGTT_i.xpt")
OGTT_all <- dplyr::bind_rows(list(OGTT_e,OGTT_f,OGTT_g,OGTT_h,OGTT_i))
OGTT_all$LBXGLT

#获取CKD数据
kiq_u_e <- read_xpt("2007-2008/Questionnaire/kiq_u_e.xpt")
kiq_u_f <- read_xpt("2009-2010/Questionnaire/kiq_u_f.xpt")
kiq_u_g <- read_xpt("2011-2012/Questionnaire/kiq_u_g.xpt")
kiq_u_h <- read_xpt("2013-2014/Questionnaire/kiq_u_h.xpt")
kiq_u_i <- read_xpt("2015-2016/Questionnaire/kiq_u_i.xpt")
kiq_u_all <- dplyr::bind_rows(list(kiq_u_e,kiq_u_f,kiq_u_g,kiq_u_h,kiq_u_i))
kiq_u_all$KIQ022
output <- plyr::join_all(list(demo_all, Laboratory_all,PHTHTE_all,BMI_all,cot_all,alq_all,alb_cr_all,bpq_all
                              ,paq_all,ghb_all,glu_all,OGTT_all,kiq_u_all),
                         by='SEQN',type='full')
#年龄筛选
output <- output[output$RIDAGEYR>=20,]
dim(output)
#尿酸数据筛选
output <- output[!is.na(output$LBDSUASI),]
dim(output)
#邻苯二甲酸酯数据筛选
output <- output[!is.na(output$URXCNP),]
output <- output[!is.na(output$URXCOP),]
output <- output[!is.na(output$URXECP),]
output <- output[!is.na(output$URXMBP),]
output <- output[!is.na(output$URXMC1),]
output <- output[!is.na(output$URXMEP),]
output <- output[!is.na(output$URXMHH),]
output <- output[!is.na(output$URXMOH),]
output <- output[!is.na(output$URXMZP),]
output <- output[!is.na(output$URXMIB),]
#BMI数据筛选
output <- output[!is.na(output$BMXBMI),]#8409
dim(output)
#家庭贫困数据筛选
output <- output[!is.na(output$INDFMPIR),]#7650
dim(output)
#饮酒数据筛选 
output <- output[!is.na(output$ALQ101),]
outnadrink1 <- which(is.na(output$ALQ101))
dim(output)
#血清可替宁数据筛选
output <- output[!is.na(output$LBXCOT),]#8407


#教育程度数据筛选
output <- output[!is.na(output$DMDEDUC2),]#7047
#尿肌酐数据筛选
output <- output[!is.na(output$URXCRS),]#7045
dim(output)
#高血压数据筛选
a <- output
output <- output[!is.na(output$BPQ020),]#7045
dim(a)
#体力活动数据筛选PAQ665
#output$sports <- ifelse(output$PAQ650==2,0,ifelse(output$PAQ665==2,1,NA))
#output <- output[!is.na(output$PAQ665),]#7045
#dim(output)
#判断糖尿病
demo1 <- output
demo1$diabetes <- ifelse(demo1$LBXGH>5.7 | demo1$LBXGLU>100 | demo1$LBXGLU>140,
                         1,0)
#性别
demo1$Sex <- ifelse(demo1$RIAGENDR==1,'male','female')
#bmi
demo1$BMI <- ifelse(demo1$BMXBMI<=24.9,'≤ 24.9', '> 24.9')

# 种族
demo1$race <- recode_factor(demo1$RIDRETH1, 
                      `1` = 'Mexican American',
                      `2` = 'Other Hispanic',
                      `3` = 'White',
                      `4` = 'Black',
                      `5` = 'Other Race'
                      
)
# 教育程度
demo1$edu <- recode_factor(demo1$DMDEDUC2, 
                     `1` = 'Less than high school',
                     `2` = 'Less than high school',
                     `3` = 'High school graduate/GED or equivalent',
                     `4` = 'College or above',
                     `5` = 'College or above'
                     
)
# 贫困程度
demo1$FPL <- ifelse(demo1$INDFMPIR<1.0,'< 1.0',
                         ifelse((demo1$INDFMPIR>=1.0&demo1$INDFMPIR<2.0) , "1.0–2.0" ,
                                ifelse((demo1$INDFMPIR>=2.0&demo1$INDFMPIR<4.0) , "2.0–4.0" ,'> 4.0')))
#活动
demo1$activity <- ifelse(demo1$PAQ605==1,'vigorous',
                              ifelse((demo1$PAQ620==1) , "moderate" ,'no'))
# 吸烟 1 是 
demo1$smoke <- ifelse(demo1$LBXCOT>10,'yes','no')
# 饮酒
demo1$drink <- ifelse(demo1$ALQ101==1,'yes','no')
#高血压
demo1$Hypertension <- ifelse(demo1$BPQ020==1,'yes','no')
# ckd 
demo1$CKD <- ifelse(demo1$KIQ022==1,'yes','no')
# 尿肌酐
demo1$urinecreatinine <- demo1$URXCRS
# 年龄
demo1$Age <- demo1$RIDAGEYR

demo1$Hyperuricemia <- ifelse((demo1$Sex=='male' & demo1$LBDSUASI>=416.0),'Yes',
                                   ifelse((demo1$Sex=='male' & demo1$LBDSUASI<416.0),'No',
                                          ifelse((demo1$Sex=='female' & demo1$LBDSUASI>=357.0),'Yes',
                                                 ifelse((demo1$Sex=='female' & demo1$LBDSUASI<357.0),'No',""))))

demo <-demo1[,c('SEQN', 'RIDAGEYR', 'Sex','Age', 'BMI','race','edu','FPL', 'activity','smoke','drink','Hypertension','diabetes','CKD',
                              'urinecreatinine','URXECP','URXMBP','URXMHH','URXMOH','URXMIB','URXCNP','URXCOP','URXMC1','URXMEP','URXMZP','LBDSUASI','Hyperuricemia')]


##################### 特征表table1 #######################
#参与的变量
myVars <- c('Sex','Age','BMI','race','edu','FPL','activity','smoke','drink','Hypertension','diabetes','CKD','urinecreatinine')
#因子变量
catVars <- c('Sex','BMI','race','edu','FPL','activity','drink','Hypertension','ckd')
#正态分布为均数+标准差,非正态分布为中位数(四分位数)
#以下为非正态分布变量
nonvar <- c('urinecreatinine')

tab <- CreateTableOne(vars = myVars, strata =c("Hyperuricemia")  , data = demo, factorVars = catVars)

tabMat <- print(tab, nonnormal = nonvar, exact = catVars,test = TRUE, smd = TRUE,showAllLevels = TRUE)
tabMat2 <- tabMat[,c(-5,-6)]
tabMat2


#######################邻苯二甲酸酯暴露水平的描述 Table S3 ####################
demoS2 <-demo[,c('URXECP', 'URXMBP', 'URXMHH', 'URXMOH','URXMIB','URXCNP','URXCOP', 'URXMC1','URXMEP','URXMZP','Hyperuricemia')]
colnames(demoS2) <- c('MECPP','MnBP','MEHHP','MEOHP','MiBP','cx-MiNP','MCOP','MCPP','MEP','MBzP','Hyperuricemia')
#参与的变量
myVars1 <- c('MECPP','MnBP','MEHHP','MEOHP','MiBP','cx-MiNP','MCOP','MCPP','MEP','MBzP')
nonvar1 <- c('MECPP','MnBP','MEHHP','MEOHP','MiBP','cx-MiNP','MCOP','MCPP','MEP','MBzP')
tab1 <- CreateTableOne(vars = myVars1, strata =c("Hyperuricemia")  , data = demoS2)
tabMat1 <- print(tab1, nonnormal = nonvar1,test = TRUE)
tabMat1 <- tabMat1[,c(-4)]
tabMat1
tab1
summary(tabMat1)


############################### 两个模型的对比 Table S4 ########################
#暂不清楚它的quantile列和PFDR列是怎么添加的,其他列都可以得到
demo$MECPP_quantile <- cut(demo$URXECP,breaks = quantile(demo$URXECP),labels = c('Q1', 'Q2', 'Q3', 'Q4'))
demo$MnBP_quantile <- cut(demo$URXMBP,breaks = quantile(demo$URXMBP),labels = c('Q1', 'Q2', 'Q3', 'Q4'))  
demo$MEHHP_quantile <- cut(demo$URXMHH,breaks = quantile(demo$URXMHH),labels = c('Q1', 'Q2', 'Q3', 'Q4'))  
demo$MEOHP_quantile <- cut(demo$URXMOH,breaks = quantile(demo$URXMOH),labels = c('Q1', 'Q2', 'Q3', 'Q4'))  
demo$MiBP_quantile <- cut(demo$URXMIB,breaks = quantile(demo$URXMIB),labels = c('Q1', 'Q2', 'Q3', 'Q4'))  
demo$cxMiNP_quantile <- cut(demo$URXCNP,breaks = quantile(demo$URXCNP),labels = c('Q1', 'Q2', 'Q3', 'Q4'))  
demo$MCOP_quantile <- cut(demo$URXCOP,breaks = quantile(demo$URXCOP),labels = c('Q1', 'Q2', 'Q3', 'Q4'))  
demo$MCPP_quantile <- cut(demo$URXMC1,breaks = quantile(demo$URXMC1),labels = c('Q1', 'Q2', 'Q3', 'Q4'))  
demo$MEP_quantile <- cut(demo$URXMEP,breaks = quantile(demo$URXMEP),labels = c('Q1', 'Q2', 'Q3', 'Q4'))  
demo$MBzP_quantile <- cut(demo$URXMZP,breaks = quantile(demo$URXMZP),labels = c('Q1', 'Q2', 'Q3', 'Q4')) 

varis<- c('MECPP_quantile','MnBP_quantile','MEHHP_quantile',
          'MEOHP_quantile','MiBP_quantile','cxMiNP_quantile','MCOP_quantile',
          'MCPP_quantile','MEP_quantile','MBzP_quantile','urinecreatinine')
demo$Hyperuricemiaint <- ifelse(demo$Hyperuricemia=='Yes',1,0)

fml<- as.formula(paste0('Hyperuricemiaint==1~',paste0(varis,collapse = '+'))) #P<0.05也是可以的
demo <- na.omit(demo)
modelA<-glm(fml,data = demo,family=binomial(link = "logit"))

model<-step(modelA,direction = "both")
glm<-summary(modelA)

OR<-round(exp(glm$coefficients[,1]),2)
SE<-round(glm$coefficients[,2],3)
CI2.5<-round(exp(coef(model)-1.96*SE),2)
CI97.5<-round(exp(coef(model)+1.96*SE),2)
CI<-paste0(CI2.5,'-',CI97.5)
B<-round(glm$coefficients[,1],3)
Z<-round(glm$coefficients[,3],3)
P<-round(glm$coefficients[,4],3)
PFDR <- p.adjust(P)

mlogit<-data.frame( 
  'B'=B,
  'SE'=SE,
  'OR'=OR,
  'CI'=CI,
  'Z' =Z,
  'P'=P)[-1,]   #-1是指删除常数项

varis1<- c('Age','Sex','race','MECPP_quantile','edu','FPL','activity','smoke',
           'drink','Hypertension','CKD','urinecreatinine','BMI','MnBP_quantile',
           'MEHHP_quantile','MEOHP_quantile','MiBP_quantile','cxMiNP_quantile',
           'MCOP_quantile','MCPP_quantile','MEP_quantile','MBzP_quantile')
fml<- as.formula(paste0('Hyperuricemiaint==1~',paste0(varis1,collapse = '+'))) #P<0.05也是可以的
demo <- na.omit(demo)
modelB<-glm(fml,data = demo,family=binomial(link = "logit"))
glm<-summary(modelB)

OR<-round(exp(glm$coefficients[,1]),2)
SE<-round(glm$coefficients[,2],3)
CI2.5<-round(exp(coef(model)-1.96*SE),2)
CI97.5<-round(exp(coef(model)+1.96*SE),2)
CI<-paste0(CI2.5,'-',CI97.5)
B<-round(glm$coefficients[,1],3)
Z<-round(glm$coefficients[,3],3)
P<-round(glm$coefficients[,4],3)
PFDR <- p.adjust(P)

mlogit1<-data.frame( 
  'B'=B,
  'SE'=SE,
  'OR'=OR,
  'CI'=CI,
  'Z' =Z,
  'P'=P)[-1,]   #-1是指删除常数项
mlogit1

######################## 邻苯二甲酸酯相关性热图 图2 ############################
#应该是数据对不上所以相关性与论文不同


demo2 <-demo1[,c('URXECP', 'URXMBP', 'URXMHH', 'URXMOH','URXMIB','URXCNP','URXCOP', 'URXMC1','URXMEP','URXMZP')]
colnames(demo2) <- c('MECPP','MnBP','MEHHP','MEOHP','MiBP','cx-MiNP','MCOP','MCPP','MEP','MBzP')

data<-cor(demo2,method="spearman")

corrplot(data, method = "number", type = "upper",order = 'hclust',
         tl.col = "black", tl.cex = 0.8, tl.srt = 45,tl.pos = "lt")
corrplot(data, method = "circle", type = "lower",order = 'hclust',
         tl.col = "n", tl.cex = 0.8, tl.pos = "n",
         add = T)
################################## MI补充数据 ###################################


demoMI <- mice(demo, #数据集
               method = "rf", #采用随机森林插补
               m=5, # 5次插补
               printFlag = FALSE #不显示历史记录
               )
demoMI <- complete(demoMI, action = 2)


##############  多变量 logistic 模型的单独敏感性分析 表S5 ######################
colnames(demoMI)
demoMI$Hyperuricemiaint <- ifelse(demoMI$Hyperuricemia=='Yes',1,0)
#补全后的数据模型
demoMIA <-demoMI[,c('URXECP', 'URXMBP', 'URXMHH', 'URXMOH','URXMIB','URXCNP','URXCOP', 'URXMC1','URXMEP','URXMZP','Hyperuricemiaint')]
colnames(demoMIA) <- c('MECPP','MnBP','MEHHP','MEOHP','MiBP','MiNP','MCOP','MCPP','MEP','MBzP','Hyperuricemiaint')
varis<- c('MECPP','MnBP','MEHHP','MEOHP','MiBP','MiNP','MCOP','MCPP','MEP','MBzP')
fml<- as.formula(paste0('Hyperuricemiaint==1~',paste0(varis,collapse = '+'))) #P<0.05也是可以的
modelA<-glm(fml,data = demoMIA,family=binomial(link = "logit"))
glm<-summary(modelA)
OR<-round(exp(glm$coefficients[,1]),2)
SE<-round(glm$coefficients[,2],3)
CI2.5<-round(exp(coef(modelA)-1.96*SE),2)
CI97.5<-round(exp(coef(modelA)+1.96*SE),2)
CI<-paste0(CI2.5,'-',CI97.5)
P<-round(glm$coefficients[,4],3)
mlogitA<-data.frame( 
  'model'='model1',
  'OR'=OR,
  'CI'=CI,
  'P'=P)[-1,]   #-1是指删除常数项

varis1<- c('Age','Sex','race','edu','FPL','activity','smoke',
           'drink','Hypertension','CKD','urinecreatinine','BMI','MECPP','MnBP','MEHHP',
           'MEOHP','MiBP','MiNP','MCOP','MCPP','MEP','MBzP')
demoMIB <-demoMI[,c('Age','Sex','race','edu','FPL','activity','smoke',
                   'drink','Hypertension','CKD','urinecreatinine','BMI','URXECP', 'URXMBP', 
                   'URXMHH', 'URXMOH','URXMIB','URXCNP','URXCOP', 'URXMC1','URXMEP','URXMZP',
                   'Hyperuricemiaint')]
colnames(demoMIB) <- c('Age','Sex','race','edu','FPL','activity','smoke',
                       'drink','Hypertension','CKD','urinecreatinine','BMI','MECPP','MnBP','MEHHP',
                       'MEOHP','MiBP','MiNP','MCOP','MCPP','MEP','MBzP','Hyperuricemiaint')
fml<- as.formula(paste0('Hyperuricemiaint==1~',paste0(varis1,collapse = '+'))) #P<0.05也是可以的

modelB<-glm(fml,data = demoMIB,family=binomial(link = "logit"))
glm<-summary(modelB)
OR<-round(exp(glm$coefficients[,1]),2)
SE<-round(glm$coefficients[,2],3)
CI2.5<-round(exp(coef(modelB)-1.96*SE),2)
CI97.5<-round(exp(coef(modelB)+1.96*SE),2)
CI<-paste0(CI2.5,'-',CI97.5)
P<-round(glm$coefficients[,4],3)
mlogitB<-data.frame( 
  'model'='model2',
  'OR'=OR,
  'CI'=CI,
  'P'=P)[-1,]   #-1是指删除常数项
mlogitB <- t(mlogitB)
mlogitB <- mlogitB[,c('MECPP','MnBP','MEHHP','MEOHP','MiBP','MiNP','MCOP','MCPP','MEP','MBzP')]
mlogitB <- t(mlogitB)
mlogitB <- as.data.frame(mlogitB)
mlogitA$x <- c('MECPP','MnBP','MEHHP','MEOHP','MiBP','MiNP','MCOP','MCPP','MEP','MBzP')
mlogitB$x <- c('MECPP','MnBP','MEHHP','MEOHP','MiBP','MiNP','MCOP','MCPP','MEP','MBzP')
MI <- rbind(mlogitA[1,],mlogitB[1,])
for (r in 2:10){
  print(r)
  print(rbind(mlogitA[r,],mlogitB[r,]))
  MI <- rbind(MI,rbind(mlogitA[r,],mlogitB[r,]))
}
MI$x <- c('MECPP','','MnBP','','MEHHP','','MEOHP','','MiBP','','MiNP','',
                 'MCOP','','MCPP','','MEP','','MBzP','')

#未补全数据模型
demo$Hyperuricemiaint <- ifelse(demo$Hyperuricemia=='Yes',1,0)
demoA <-demo[,c('URXECP', 'URXMBP', 'URXMHH', 'URXMOH','URXMIB','URXCNP','URXCOP', 'URXMC1','URXMEP','URXMZP','Hyperuricemiaint')]
colnames(demoA) <- c('MECPP','MnBP','MEHHP','MEOHP','MiBP','MiNP','MCOP','MCPP','MEP','MBzP','Hyperuricemiaint')
varis<- c('MECPP','MnBP','MEHHP','MEOHP','MiBP','MiNP','MCOP','MCPP','MEP','MBzP')
fml<- as.formula(paste0('Hyperuricemiaint==1~',paste0(varis,collapse = '+'))) #P<0.05也是可以的
modelC<-glm(fml,data = demoA,family=binomial(link = "logit"))
glm<-summary(modelC)
OR<-round(exp(glm$coefficients[,1]),2)
SE<-round(glm$coefficients[,2],3)
CI2.5<-round(exp(coef(modelC)-1.96*SE),2)
CI97.5<-round(exp(coef(modelC)+1.96*SE),2)
CI<-paste0(CI2.5,'-',CI97.5)
P<-round(glm$coefficients[,4],3)
mlogitC<-data.frame( 
  'model'='model1',
  'OR'=OR,
  'CI'=CI,
  'P'=P)[-1,]   #-1是指删除常数项

varis1<- c('Age','Sex','race','edu','FPL','activity','smoke',
           'drink','Hypertension','CKD','urinecreatinine','BMI','MECPP','MnBP','MEHHP',
           'MEOHP','MiBP','MiNP','MCOP','MCPP','MEP','MBzP')
demoB <-demo[,c('Age','Sex','race','edu','FPL','activity','smoke',
                'drink','Hypertension','CKD','urinecreatinine','BMI','URXECP', 'URXMBP', 
                'URXMHH', 'URXMOH','URXMIB','URXCNP','URXCOP', 'URXMC1','URXMEP','URXMZP',
                'Hyperuricemiaint')]
colnames(demoB) <- c('Age','Sex','race','edu','FPL','activity','smoke',
                     'drink','Hypertension','CKD','urinecreatinine','BMI','MECPP','MnBP','MEHHP',
                     'MEOHP','MiBP','MiNP','MCOP','MCPP','MEP','MBzP','Hyperuricemiaint')
fml<- as.formula(paste0('Hyperuricemiaint==1~',paste0(varis1,collapse = '+'))) #P<0.05也是可以的

modelD<-glm(fml,data = demoB,family=binomial(link = "logit"))
glm<-summary(modelD)
OR<-round(exp(glm$coefficients[,1]),2)
SE<-round(glm$coefficients[,2],3)
CI2.5<-round(exp(coef(modelD)-1.96*SE),2)
CI97.5<-round(exp(coef(modelD)+1.96*SE),2)
CI<-paste0(CI2.5,'-',CI97.5)
P<-round(glm$coefficients[,4],3)
mlogitD<-data.frame( 
  'model'='model2',
  'OR'=OR,
  'CI'=CI,
  'P'=P)[-1,]   #-1是指删除常数项
mlogitD <- t(mlogitD)
mlogitD <- mlogitD[,c('MECPP','MnBP','MEHHP','MEOHP','MiBP','MiNP','MCOP','MCPP','MEP','MBzP')]
mlogitD <- t(mlogitD)
mlogitD <- as.data.frame(mlogitD)
mlogitC$x <- c('MECPP','MnBP','MEHHP','MEOHP','MiBP','MiNP','MCOP','MCPP','MEP','MBzP')
mlogitD$x <- c('MECPP','MnBP','MEHHP','MEOHP','MiBP','MiNP','MCOP','MCPP','MEP','MBzP')
Subgroup <- rbind(mlogitC[1,],mlogitD[1,])
for (r in 2:10){
  print(r)
  print(rbind(mlogitC[r,],mlogitD[r,]))
  Subgroup <- rbind(Subgroup,rbind(mlogitC[r,],mlogitD[r,]))
}
Subgroup
Subgroup$x <- c('MECPP','','MnBP','','MEHHP','','MEOHP','','MiBP','','MiNP','',
          'MCOP','','MCPP','','MEP','','MBzP','')

cbind(MI,Subgroup)
##################### 限制性立方样条(RCS) Table S2 S3 ##########################
#暂时先这样,后续看看还有没有其他方法,或者看论文里有没有更具体的关于RCS使用参数的
#解释,目前感觉就是x轴的值对不上,应该是对变量进行了一次处理然后再进行的建模和绘制图形
#先跳过吧!!!!!!!!!!!!
demoB <-demo[,c('Age','Sex','race','edu','FPL','activity','smoke',
                'drink','Hypertension','CKD','urinecreatinine','BMI','URXECP', 'URXMBP', 
                'URXMHH', 'URXMOH','URXMIB','URXCNP','URXCOP', 'URXMC1','URXMEP','URXMZP',
                'Hyperuricemiaint','LBDSUASI')]
colnames(demoB) <- c('Age','Sex','race','edu','FPL','activity','smoke',
                     'drink','Hypertension','CKD','urinecreatinine','BMI','MECPP','MnBP','MEHHP',
                     'MEOHP','MiBP','MiNP','MCOP','MCPP','MEP','MBzP','Hyperuricemiaint','LBDSUASI')

#MINP  S2
ddist<-datadist(demoB)
options(datadist="ddist")
fit <-lrm(Hyperuricemiaint ~rcs(MiNP,nk=5)+Age + Sex + race + edu + FPL + activity + 
            smoke + drink + Hypertension + CKD + urinecreatinine + BMI + 
            MECPP + MnBP + MEHHP + MEOHP + MiBP + MCOP + MCPP + 
            MEP + MBzP,data=demoB,x=FALSE)
pred_OR<-Predict(fit,MiNP,ref.zero=TRUE,fun=exp)
par(mar = c(5, 4, 4, 4) + 0.3)
par(xpd=NA)
ylim.bot<-min(pred_OR[,"lower"])
ylim.top<-max(pred_OR[,"upper"])
plot(pred_OR[,1],pred_OR[,"yhat"], 
     xlab = "MiNP",ylab = paste0("OR"),
     type = "l",ylim = c(ylim.bot,ylim.top),
     col="red",lwd=2) 
lines(pred_OR[,1],pred_OR[,"lower"],lty=2,lwd=1.5)
lines(pred_OR[,1],pred_OR[,"upper"],lty=2,lwd=1.5)
lines(x=range(pred_OR[,1]),y=c(1,1),lty=3,col="grey40",lwd=1.3) #
points(1,1,pch=16,cex=1.2)
#MCOP S2
ddist<-datadist(demoB)
options(datadist="ddist")
fit <-lrm(Hyperuricemiaint ~rcs(MCOP,nk=5)+Age + Sex + race + edu + FPL + activity + 
            smoke + drink + Hypertension + CKD + urinecreatinine + BMI + 
            MECPP + MnBP + MEHHP + MEOHP + MiBP + MiNP + MCPP + 
            MEP + MBzP,data=demoB,x=FALSE)
pred_OR<-Predict(fit,MCOP,ref.zero=TRUE,fun=exp)
par(mar = c(5, 4, 4, 4) + 0.3)
par(xpd=NA)
ylim.bot<-min(pred_OR[,"lower"])
ylim.top<-max(pred_OR[,"upper"])
plot(pred_OR[,1],pred_OR[,"yhat"], 
     xlab = "MCOP",ylab = paste0("OR"),
     type = "l",ylim = c(ylim.bot,ylim.top),
     col="red",lwd=2) 
lines(pred_OR[,1],pred_OR[,"lower"],lty=2,lwd=1.5)
lines(pred_OR[,1],pred_OR[,"upper"],lty=2,lwd=1.5)
lines(x=range(pred_OR[,1]),y=c(1,1),lty=3,col="grey40",lwd=1.3) #
points(2.5,1,pch=16,cex=1.2)


#MINP  S3
ddist<-datadist(demoB)
options(datadist="ddist")
demoB$LBDSUASI
fit <-lrm(LBDSUASI ~rcs(MiNP,nk=5)+Age + Sex + race + edu + FPL + activity + 
            smoke + drink + Hypertension + CKD + urinecreatinine + BMI + 
            MECPP + MnBP + MEHHP + MEOHP + MiBP + MCOP + MCPP + 
            MEP + MBzP,data=demoB,x=FALSE)
pred_OR<-Predict(fit,MiNP,ref.zero=TRUE,fun=exp)
par(mar = c(5, 4, 4, 4) + 0.3)
par(xpd=NA)
ylim.bot<-min(pred_OR[,"lower"])
ylim.top<-max(pred_OR[,"upper"])
plot(pred_OR[,1],pred_OR[,"yhat"], 
     xlab = "MiNP",ylab = paste0("OR"),
     type = "l",ylim = c(ylim.bot,ylim.top),
     col="red",lwd=2,xlim = c(-2,5)) 
lines(pred_OR[,1],pred_OR[,"lower"],lty=2,lwd=1.5)
lines(pred_OR[,1],pred_OR[,"upper"],lty=2,lwd=1.5)
lines(x=range(pred_OR[,1]),y=c(1,1),lty=3,col="grey40",lwd=1.3)
points(1,1,pch=16,cex=1.2)


########################## 贝叶斯模型及绘制图片 ###############################
#引入贝叶斯r包# BKMR公式中的x为协变量，z才是暴露因素

# res <- persp(z1, z2, hgrid.true, theta = 30, phi = 20, expand = 0.5, 
#              col = "lightblue", xlab = "", ylab = "", zlab = "")
# 小批量数据测试
demoBTest <- demoB[c(1:200),]

y <- demoBTest$Hyperuricemiaint
Z <- demoBTest[,c('MECPP','MnBP','MEHHP','MEOHP','MiBP','MiNP',
              'MCOP','MCPP','MEP','MBzP')]
#X <- demoB[,c('Age','Sex','race','edu','FPL','activity','smoke',
#              'drink','Hypertension','CKD','urinecreatinine','BMI')]
X <- demoBTest[,c('Age','urinecreatinine')]
#先使用模型一试试看,因为X必须是number
#fitkm <- kmbayes(y =y, Z = Z, X = X, iter = 10000, verbose = FALSE, varsel = TRUE)
fitkm <- kmbayes(y =y, Z = Z, X = NULL, iter = 5000,verbose = FALSE, varsel = TRUE, family="binomial",est.h=TRUE)
# str(X)

pred.resp.univar <- PredictorResponseUnivar(fit = fitkm)
# ggplot(pred.resp.univar, aes(z, est, ymin = est - 1.96*se, ymax = est + 1.96*se)) + 
ggplot(pred.resp.univar, aes(z, est, ymin = est - 1.96*se, ymax = est + 1.96*se)) + 
  geom_smooth(stat = "identity") + 
  facet_wrap(~ variable) +
  ylab("h(expos)") +
  xlab("Phthalate metabolites")

#####
pred.resp.bivar <- PredictorResponseBivar(fit = fitkm, min.plot.dist = 1)
pred.resp.bivar.levels <- PredictorResponseBivarLevels(
  pred.resp.df = pred.resp.bivar, 
  Z = Z, qs = c(0.1, 0.5, 0.9))
ggplot(pred.resp.bivar.levels, aes(z1, est)) + 
  geom_smooth(aes(col = quantile), stat = "identity") + 
  facet_grid(variable2 ~ variable1) +
  ggtitle("h(expos1 | quantiles of expos2)") +
  xlab("expos1")

#####
risks.overall <- OverallRiskSummaries(fit = fitkm, y = y, Z = Z, X = NULL, 
                                      qs = seq(0.25, 0.75, by = 0.05), 
                                      q.fixed = 0.5, method = "exact")
risks.overall
ggplot(risks.overall, aes(quantile, est, ymin = est - 1.96*sd, ymax = est + 1.96*sd)) + 
  geom_pointrange()


#####
y1 <- demoB[,'LBDSUASI']
fitkm <- kmbayes(y =y1, Z = Z, X = NULL, iter = 2,verbose = FALSE, varsel = TRUE)
risks.singvar <- SingVarRiskSummaries(fit = fitkm, y = y1, Z = Z, X = NULL, 
                                      qs.diff = c(0.25, 0.75), 
                                      q.fixed = c(0.25, 0.50, 0.75),
                                      method = "exact")
risks.singvar

ggplot(risks.singvar, aes(variable, est, ymin = est - 1.96*sd, 
                          ymax = est + 1.96*sd, col = q.fixed)) + 
  geom_pointrange(position = position_dodge(width = 0.75)) + 
  coord_flip()










